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We propose a new algorithm for non-equilibrium molecular dynamics simulations of thermal gradients. The algorithm 
is an extension of the heat exchange algorithm developed by Hafskjold and co-workers [Mol. Phys. 80 , 1389 (1993); 
Mol. Phys. 81 , 251 (1994)], in which a certain amount of heat is added to one region and removed from another by 
rescaling velocities appropriately. Since the amount of added and removed heat is the same and the dynamics between 
velocity rescaling steps is Hamiltonian, the heat exchange algorithm is expected to conserve the energy. However, it has 
been reported previously that the original version of the heat exchange algorithm exhibits a pronounced drift in the total 
energy, the exact cause of which remained hitherto unclear. Here, we show that the energy drift is due to the truncation 
error arising from the operator splitting and suggest an additional coordinate integration step as a remedy. The new 
algorithm retains all the advantages of the original one whilst exhibiting excellent energy conservation as illustrated for a 
Lennard-Jones liquid and SPC/E water. 
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I. INTRODUCTION 

Non-equilibrium molecular dynamics (NEMD) simulations 
allow us to study transport phenomena and determine transport 
coefficients'. In studies of heat conduction, an external field is 
applied to the system thereby driving it to a steady state. The 
nature of the coupling between the external field and the system 
differs between algorithms and determines whether a spatially 
homogeneous state^, a temperature gradient^^^ or a heat flux is 
imposed®^®. A suitable algorithm for a particular application 
depends on its ability to model the underlying physics correctly. 
If energy is supplied at a constant rate in an experiment, for 
example, a thermostat which imposes a heat flux would lend 
itself for the simulation. Erom a computational point of view, 
generating a flux might be preferable, because it is simpler to 
measure the temperature than the heat flux. 

One way to generate a heat flux in computer simulations 
involves swapping kinetic energy between two subdomains of 
the simulation box®’^. In the heat exchange (HEX) algorithm 
developed by Ikeshoji and Hafskjold®’*", a specific amount of 
heat is periodically removed from one subdomain or reservoir, 
and supplied to the other. These two regions thus act as a heat 
sink and source, respectively. The HEX method adjusts the 
non-translational kinetic energy by velocity rescaling while 
preserving the individual center of mass velocities of the two 
heat reservoirs. Other methods use different procedures to gen¬ 
erate heat fluxes. In the reverse NEMD (RNEMD) method 
developed by Miiller-Plathe^, the heat transfer is established 
by continuously identifying hot and cold particles inside the 
reservoirs and exchanging their momenta. Extensions of the 
RNEMD method were proposed by Kuang and Gezelter**’", 
who replaced the momenta swaps by velocity rescaling moves. 
The velocity scaling and shearing (VSS) RNEMD method" al¬ 
lows for imposing a momentum flux in addition to the thermal 
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flux. However, we note that in the absence of any momentum 
flux, the method is identical to the HEX algorithm. Although 
these methods are widely applicable, they all lack an attrac¬ 
tive feature which is a formulation based on time continuous 
equations of motion. Knowing the equations of motion is ad¬ 
vantageous, for example, if one is interested in studying system 
properties such as phase space compressibility or the develop¬ 
ment of accurate integration schemes. 

Due to its simplicity, the HEX algorithm is an attractive 
choice for simulating a fluid in the absence of solid inhomo¬ 
geneities. Since the same amount of energy is added and re¬ 
moved, one would expect the algorithm to conserve the total 
energy exactly. However, as pointed out in one of the original 
papers'" and subsequent work", numerical implementations 
of the algorithm lead to a considerable energy drift over sim¬ 
ulation time scales of a few nanoseconds. A change in total 
energy of several percent of the initial value was considered 
acceptable in past work. Nevertheless, the energy drift is a 
severe restriction limiting the accessible simulation time scales. 
Remedies to this problem either involve employing a smaller 
timestep or compensating the energy drift with an additional 
thermostat*^ which is undesirable, because such thermostats 
may affect the very temperature profile that one aims to study. 

In this work, we identify the underlying cause of the energy 
loss and suggest a new algorithm to achieve improved energy 
conservation. The paper is organised as follows: In Sec. II, we 
summarise the HEX algorithm and its numerical implemen¬ 
tation. We then show that the integration scheme leads to an 
unphysical energy drift. In order to understand the origin of 
this problem, we derive the equations of motion for continu¬ 
ous time in Sec. III. This allows us to express the integration 
scheme as a Trotter factorisation of the Liouville operator. In 
Sec. IV, we work out the leading-order error term of the em¬ 
ployed operator splitting. Based on our analysis, we propose 
an enhanced algorithm in Sec. V and compare the results in 
Sec. VI. 
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where the rescaling factor is given by 
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Here, updated quantities are denoted with an overbar. It can 
easily be verified that the update step given by Eq. (2) satisfies 
iCpt. — • Since there is no net 

energy flux into the system according to our assumptions, this 
also implies that the total system energy E remains constant. 

We note that the above formulation of the velocity update as 
presented by Aubry et al. is simpler than the one which was 
originally proposed by Ikeshoji and Hafskjold®. In the latter 
case, is a more complex function of the velocities, but it is 
easy to see that both formulations are equivalent. 


FIG. 1. Illustration of the simulation box, 1?, with Hamiltonian re¬ 
gions, /o, a hot region, A (red), and a cold region, /2 (blue). The 
centre of mass velocities of 17, A and T 2 are vn, va and vr^, re¬ 
spectively. Atoms are represented by red/blue circles, if they are 
located in the hot/cold region and by empty circles otherwise. 


II. HEX ALGORITHM 

The goal of the HEX algorithm is to impose a constant heat 
flux onto the system. This is accomplished by adding heat 
AQa at each timestep to Np pair-wise disjoint subdomains 
Ek, of the simulation box 17 (Fig. 1). Heat is subtracted if 
AQ Pf, is negative. We label those parts of the simulation box 
which are not thermostatted with Jq. The box contains N 
atoms each labelled with a unique index. If there is no net 
energy flux into the simulation box as we assume here, i.e. 

AQa = 0’ the system will approach a steady state in 
which heat fluxes are established between the subdomains. The 
position and velocity vectors of atom i are and Vi, respec¬ 
tively. Furthermore, we use vp^. and ttr? to denote the centre of 
mass velocities of the regions /jj and the box 17, respectively. 


A. Energy supply 


Energy is added or removed by rescaling the velocities of 
all particles contained in region E}- by the same factor and 
shifting them by a constant. The value of is chosen such 
that the non-translational kinetic energy of that region. 
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changes by AQp^ leaving vp,, unchanged, where mp^ is the 
total mass contained in /\. The time-dependent index set 
7 fc comprises all particles which are located in 7^. Particles 
outside any thermostatted region are not affected by this proce¬ 
dure. For the individual region Ek, the velocity update can be 
formulated as®-'^ 


( 2 ) 


B. Time integration 

In order to keep track of the time evolution, it is convenient 
to introduce some additional notation. We label all quantities 
sampled at time t = nAt with a superscript n, where At is 
the timestep. In addition we define ^0 to be unity at all times 
and k{ri) to be the index of the region in which particle i is 
located. The current state of the system is fully described by 
a 6A-dimensionaI vector x = (r,r;) in phase space, where 
the vectors r and v contain all particle positions and velocities, 
respectively. The HEX algorithm for velocity Verlet^^ can then 
be formulated as 
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where U{r) is the potential energy and fi the force acting on 
particle i. For the entire scheme to be symmetric, half the en¬ 
ergy is supplied at the beginning of the timestep and the other 
half at the end. The scaling factors and ^re eval¬ 
uated using Eq. (3) at the states and 

respectively. 

We note that in the original work^’'°, the authors do not pro¬ 
vide any details about when exactly the thermostatting step 
should happen. For comparison, we also tested an asymmetric 
version of the algorithm (HEX/a), where all the energy is sup¬ 
plied at the end of the timestep. In this case, the initial velocity 
update reduces to the identity operation. 


C. Model system 

We studied the energy conservation of the HEX algorithm 
for a Lennard-Jones (LJ) fluid using the simulation package 


A = C/cA + (1 - ^k)vp^, 




















3 


LAMMPS (version 906014)'“*. The symmetric, pairwise LJ 
potential is given by*^ 


uu{r) = 4e 



where e is the depth of the potential and a the effective atomic 
diameter. In order to rule out any effects due to simple spherical 
truncation of the potential, we employed a slightly modified 
potential which is given by*® 

MsF(f’) = uuir) - ulj(a) - {r- r^)uij{r^) (6) 

forr < Ts and zero otherwise, where Tj is the cutoff. From 
the functional form of Eq. (6) it is clear that uspir) and u'^p{r) 
are both continuous at the cutoff. We employed a value of 
r* = 3 for all simulations in this section. (Reduced quantities 
are labelled with an asterisk.) 


D. Equilibration 



AiVlO”^ 

FIG. 2. Energy loss for LJ at the final time t* = 5000 as a function 
of the timestep. Each point in the figure corresponds to a separate 
simulation. The equilibrium run (blue diamonds) is compared to the 
symmetric (red circles) and asymmetric (black squares) versions of 
the HEX algorithm, respectively. A quadratic fit (red, solid line) was 
carried out for the symmetric version. 


The rectangular simulation box with dimensions T*/2 = 
L* — L* = 10.58 comprised N = 2000 atoms resulting in a 
density of p* = 0.8444. The thermodynamic conditions con¬ 
sidered in this work are similar to those in Ref. 6. Starting from 
an initial lattice structure with zero linear momentum, the sys¬ 
tem was heated up to twice the target temperature of T* = 0.72 
and subsequently cooled down again at the same rate. The 
thermo staking during this initial period was achieved by veloc¬ 
ity rescaling and the entire annealing process took 2.5 x 10^ 
timesteps. The equations of motion were integrated with the 
velocity Verlet algorithm using a timestep of At* = 0.002. We 
then increased the timestep to At* = 0.004 and carried out a 
2 X 10® timestep NVT simulation using a Nose-Hoover ther¬ 
mostat'^** with a relaxation time of r* = 0.5. During this 
run we computed the average system energy. Using the HEX 
algorithm, we then adjusted the energy of the last configura¬ 
tion and used it as input for another 2 x 10® timestep NVE 
equilibration run. This procedure allowed us to achieve an 
average equilibrium temperature of T* = (0.7200 ± 0.0002). 
The error bar corresponds to one standard deviation of the error 
of the mean, the variance of which was estimated using block 
average analysis'®. 

As a reference for the energy conservation in equilibrium, 
we carried out an additional set of NVE simulations at vari¬ 
ous timesteps. With the above protocol we matched the tem¬ 
perature of these runs to be close to the one inside the hot 
reservoir in the NEMD case. The average temperature was 
T* = (0.8400 ±0.0002). 

E. Energy conservation during NEMD 

The previously equilibrated structures were subjected to a 
temperature gradient along the z-axis using the HEX algo¬ 
rithm. Always starting from the same phase space point, we 
varied the timestep for a fixed energy flux J^Pk = A / At 

into the reservoir. The two thermostatted regions are centred 
at the points z = and have a width of 2 in reduced 


units (Fig. 1). During each timestep, the heat AQ is taken 
from ±2 and added to Ei (AQp^ = —AQp^ = AQ > 0). We 
waited for 100 reduced time units for any transient behaviour to 
disappear and to allow the system to reach a steady state. The 
production run of 5000 reduced time units started at t* = 0. 

In order to capture the spatial variation of the temperature, 
we divided the z-axis into N\, bins. We use the notation Xj for 
the evaluation of a quantity X over bin j and assign the value 
to the centre of the bin. The instantaneous kinetic temperature 
of bin j is then given by 

2JC, 

where ICj is the total non-translational kinetic energy of the 
bin, Nj the number of atoms contained in the bin and ks Boltz¬ 
mann’s constant. The quantity / is the number of degrees of 
freedom per atom (/lj = 3 and /spc/e = 2). We subtracted 
three degrees of freedom to account for the centre of mass ve¬ 
locity of the bin. In the stationary state, the heat flux between 
the reservoirs in Fig. 1 is given by 

^ AQ ^ 

2AtL^Ly 2L^Ly' ’ 

where the factor of 2 in the denominator accounts for the pe¬ 
riodic setup. Considering a reference layer, this is intuitively 
clear, because half the supplied heat will flow to the other reser¬ 
voir in the reference box and the other half to its image in the 
neighbouring box. The heat flux is an input parameter of the 
HEX algorithm, which we set to 0.15 in reduced units. 

The dependence of the energy loss at the final time on the 
timestep is shown in Eig. 2. Erom the quadratic fit, it is clear 
that the HEX algorithm exhibits an energy drift which scales 
as 0{At^). On the other hand, the energy was conserved per¬ 
fectly well in NVE simulations at the peak temperature inside 
the hot reservoir. (The temperature profiles are discussed in 
Sec. VI.) 
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III. EQUATIONS OF MOTION 


IV. OPERATOR SPLITTING 


To gain a better understanding of the energy drift of the HEX 
algorithm, we first derive the ordinary differential equations 
(ODEs) solved by the algorithm in the limit At —> 0. To 
this end we consider the velocity update for continuous time. 
Dropping all particle and region indices for readability and 
eliminating the intermediate velocities, we can cast Eq. (4f) 
into 
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If we subtract t;" on both sides and divide by the timestep, we 
get 
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Our goal is to show that the energy drift is caused by higher- 
order truncation terms, which are not taken into account in the 
time integration. These terms can be derived easily once the in¬ 
tegration scheme is expressed in terms of a Trotter factorisation 
of the Liouville operator. 


A. Trotter factorisation 

Tuckerman et al .showed that reversible integrators can 
be generated based on a Trotter factorisation of the Liouville 
operator iL. Utilising the same theoretical framework, we 
consider the splitting 


iL = iLi + iL2, 
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It is straightforward to show that 

(e+i^ - 1) v” _ Trv{t^) 

Af ^ 2/Cr(f") 


and apply it to the current state of the system which is fully 
described by x in the 6A-dimensional phase space. The exact 
time evolution of the system is formally given by 

a5ex(f) = e**-^x(0). (14) 


and 


^„+i + ^"+1) ^^+1 ^ Frvr{t^) 

At ^ 2/Cr(<") 

in the limit of Af —> 0. From Eq. (4c), it is immediately 
obvious that the derivative of the coordinates is given by the 
velocities. 

The continuous equations of motion solved by the HEX 
algorithm are therefore given by 

ri = Vi, (11a) 

• f i I /111 \ 

v, = — H-, (lib) 

nii nii 

where the thermostatting force is defined as 


Unfortunately, it is not feasible to evaluate this expression 
analytically for the equations we are interested in. The problem 
can be simplified, however, by considering the approximation 


x{t) 
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where P is an integral number which implicitly defines the 
timestep through At = t/P. The operator eT“®^i acts on 
the velocities and adds the energy AQI2 to the system. In 
fact, as shown in Appendix A, the velocity update of the HEX 
algorithm is the exact solution of this operation. Hamilton’s 
equations of motion are then integrated with followed 

by the second energy supply. For the analysis in the next 
subsection, we assume that all operations in Eq. (15) can be 
carried out analytically, although in the simulation we use an 
additional approximation of discussed in Sec. V. 
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if k{ri) > 0, 
otherwise. 


( 12 ) 


In order for the equations to be well-defined, we assume that 
there are sufficiently many particles inside any reservoir, i.e. 
regions with k{ri) > 0, such that the non-translational kinetic 
energy never vanishes. Outside the reservoirs the thermostat¬ 
ting force is zero and the particles obey Hamiltonian motion. 
Some further properties of the equations are analysed in Ap¬ 
pendix A. 


B. Local truncation error 

The splitting given by Eq. (15) is known as Strang splitting^®. 
It has the local truncation error^' 

x{At) — atex(Af) = At^£x^^{Q) + O (Af^), (16) 

where the first term on the RHS is determined by the operator 

£ = ^[iL2,[iL2,iLi]] - ^ [iLi, [zLi, iLa] ] (17) 
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and [A, B] = AB — BA is the commutator. Rearranging terms, A. Rigid molecules 
we find 


x{At) — At^£xe^{0) = Xex{At) + O (At^) . (18) 

This means that the key to improving the accuracy of the nu¬ 
merical approximation is to apply the correction —Af^5a;ex(0) 
to the original solution. Alternatively, we can also use a correc¬ 
tion —At^£x{At), where x{At) = Xex(O) + 0{At), without 
changing the order of the truncation error. 


V. ENHANCED HEX ALGORITHM 

The analysis of the previous section remains valid for any 
approximation of which is sufficiently accurate. This 

is necessarily the case if the local truncation error is O (Af^) 
or higher. Velocity Verlet integration is less accurate than that 
and has a local truncation error of O (Af^). Nevertheless, we 
found that it is fully sufficient to consider a coordinate correc¬ 
tion of the form of Eq. (18) to get hold of the energy loss. We 
therefore ignored the additional velocity Verlet tmncation error 
and all other correction terms in Eq. (18) affecting velocities 
only. 

This analysis leads us directly to the enhanced heat exchange 
(eHEX) algorithm, which is defined through the update se¬ 
quence 
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Apart from the last integration step and some relabelling, this 
scheme is identical to the HEX algorithm. As shown in Ap¬ 
pendix A, the correction term is given by 


f' T • 

^' i,oc 
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and evaluated at the state We note that this expression 

vanishes for particles outside any reservoir, because the ther¬ 
mostatting force is zero in that case. The scaling factors ) 

and calculated at the system states x^ and re¬ 

spectively. As in the formulation of the original algorithm, we 
also consider the case where all the energy is supplied asym¬ 
metrically in Eq. (19f). We refer to this version of the algorithm 
as eHEX/a. 


Employing constraining forces, we can extend the eHEX 
algorithm to a system of rigid bodies, such as SPC/E water^^. 
In the SHAKE algorithm^^, originally devised for Verlet in¬ 
tegration, rigidity is imposed by solving iteratively for a set 
of Lagrange multipliers. If the underlying equations are in¬ 
tegrated with the velocity Verlet algorithm, a second set of 
constraining forces is required to eliminate velocity compo¬ 
nents along any fixed bond. This is taken into account by the 
RATTLE algorithm^"^ which we implemented in LAMMPS. 

To be compatible with the treatment of constraining forces 
in LAMMPS, we consider the eHEX/a algorithm for rigid bod¬ 
ies. We use RATTLE to ensure that the velocities and positions 
are satisfied up to the target tolerance after the second veloc¬ 
ity update (Eq. (19e)). Provided that all sites of a reference 
molecule are located in the same region, the scaling and shift¬ 
ing in Eq. (19f) does not violate the constraints. For this reason, 
we only rescale an individual site of a molecule if its centre of 
mass is located within the reservoir. 

For the small fraction of molecules inside a reservoir, the 
coordinate correction in Eq. (19g) introduces an 0{At^) error 
in the bond distances. This error is small and of the same order 
as the local error of the RATTLE algorithm itself^^. For this 
reason, we consider an unconstrained update acceptable. How¬ 
ever, we monitored the maximum relative errors throughout all 
simulations. The constraining forces for the coordinates are re¬ 
calculated at the end of the timestep to ensure that the positions 
are correct after the subsequent velocity Verlet update. 


B. Model system 

In addition to the monatomic system, we tested the eHEX/a 
algorithm for the SPC/E water model. The simulation box 
with dimensions = Ly = 25.26 A contained 1024 

molecules resulting in a density of 0.95 g/cm^. We used a 
real-space cutoff of 11 A for the LJ and Coulomb interactions, 
which were evaluated with standard Ewald summation'®. The 
damping parameter was a = 6.816/La; with 9841 fe-vectors 
(before employing symmetry properties of the reciprocal sum). 

Starting from a lattice structure, we rescaled velocities for 
10 ps to drive the system close to a target temperature of 400 K. 
We employed a timestep of 1 fs and the equations were inte¬ 
grated with velocity Verlet. This was followed by a 2 x 10® 
timestep NVT simulation using a Nose-Hoover thermostat with 
a relaxation time of 1 ps. The total energy of the last configura¬ 
tion was then adjusted such that it corresponds to the average 
of the NVT run. The average temperature over a subsequent 
2 X 10® timestep WE simulation was (400.5 ± 0.2) K. 

We then switched on the thermostat and waited for 100 ps 
for the system to reach a steady state before starting with the 
1 ns production run. The reservoirs of width 4 A were centred 
at the points z = and we imposed a heat flux of 4.08 x 

10'" W/m^. As a reference for the energy conservation, we 
carried out an additional set of 1 ns NVE simulations at various 
timesteps. The temperature in these mns was (468.0 ± 0.2) K. 
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VI. RESULTS 

The effect of the additional coordinate integration in the 
eHEX algorithm on the total energy conservation is shown in 
Figs. 3-4. As can be seen, the new algorithm exhibits excellent 
energy conservation. Even for large timesteps close to the 
stability limit for a NVE simulation at the peak temperature 
LJ ~ 0.0075 and Afmax.spc/E « 3.5 fs), there is no 
noticeable drift on this scale. The energy loss of the HEX 
algorithm, on the other hand, is substantial. At the largest 
timestep, the total system energy changed by about 0.45% for 
LJ and 1.6% for SPC/E, respectively. Although an energy loss 
of several percent was considered acceptable in the past®’'°’'\ it 
sets an upper limit to the accessible simulation time scales. The 
only way to circumvent this problem apart from coupling the 
system to an additional thermostat is to decrease the timestep 
and thereby waste valuable computing time. Based on a series 
of eight simulations at the largest timestep and with different 
initial conditions, we can give a conservative estimate of the 
improvement due to the new algorithms. For LJ we found that 
the eHEX algorithm loses at least 500 times less energy than 
the HEX algorithm (450 for eHEX/a as compared to HEX/a). 
For SPC/E water we found that the eHEX/a algorithm improves 
the energy conservation by at least a factor of 100 as compared 
to the HEX/a algorithm. The accessible simulation time scale 
therefore increased by two orders of magnitude. 

The spatial variation in temperature is shown in Figs. 5-6. 
For the monatomic system, the results agree well without any 
marked differences (Fig. 5). We note that in both cases there 
are some visible discontinuities in the vicinity of the reservoirs. 
This should not be very surprising, since the thermostatting 
force is also discontinuous. We found that the gap decreases 
as we go to lower temperature gradients, because the fluid can 
dissipate the heat sufficiently fast. A possible way of control¬ 
ling this gap is to employ a weight function and to redefine 
temperature such that it is continuous at the boundary of the 
reservoir'*. This procedure allows for better control and is 
numerically convenient, but it is not obvious which weight 
function is physically most meaningful. Furthermore, general¬ 
ising this approach to arbitrary reservoir shapes is challenging, 
because it requires some sort of signed distance information to 
the boundary. 

For SPC/E water the energy loss at large timesteps is re¬ 
flected in a slight drop of temperature (Fig. 6). The overall 
profiles agree well, but they are shifted by a few Kelvin. This 
shift is consistent with the energy loss of about 1% for the 2.5 fs 
timestep. There are no visible temperature discontinuities in 
the vicinity of the reservoirs. This might be related to the fact 
that in our scheme the boundaries are naturally smeared out as 
we only rescale entire molecules which could be intersected by 
the reservoir boundary. 

Although we omitted a constrained coordinate update in the 
eHEX/a algorithm, the relative deviation from the ideal bond 
distance never exceeded 1.1 x 10“^. This was the case for the 
largest timestep of 3 fs, but the error decays rapidly (with Af^) 
such that it reduced to 3.6 x 10“*^ for a timestep of 1 fs. The 
maximum induced relative velocity along any rigid bond was 
an order of magnitude lower for both timesteps, respectively. 
Only a small fraction of molecules inside a reservoir (« 16%) 



AC/IO”^ 

FIG. 3. Energy loss for LJ at the final time t* = 5000 for various 
timesteps. The equilibrium run (blue diamonds) and the symmetric 
(black squares) and asymmetric (green, open circles) versions of the 
eHEX algorithm, respectively, do not show any appreciable drift. The 
energy loss of the HEX algorithm (red, full circles) together with a 
quadratic fit (red, solid line) is shown for comparison. 



Ai/fs 

EIG. 4. Energy loss for SPC/E water at the final time f = 1 ns for 
various timesteps. The equilibrium run (blue diamonds) is compared 
to the asymmetric eHEX algorithm (green, open circles) and the asym¬ 
metric HEX algorithm (red, full circles) together with a quadratic fit 
(red, solid line). 


suffers from this inconsistency. We consider this error accept¬ 
able and an unconstrained update justified. An extension of 
the eHEX algorithm to a constrained update is possible in case 
higher precision is required. 

With regard to conservation of total linear momentum, we 
found that both algorithms satisfied this condition perfectly. 
We initialised the linear momentum of the box to zero at the be¬ 
ginning and it remained close to machine precision throughout 
the entire simulation. 


VII. CONCLUSIONS 

In this paper, we have presented a new algorithm for NEMD 
simulations of thermal gradients. The method comprises an 
extension to the HEX algorithm, which rescales and shifts 
velocities of particles inside reservoirs to impose a constant 
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FIG. 5. Comparison of the temperature profiles for LJ. 



FIG. 6. Comparison of the temperature profiles for SPC/E water. 
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Appendix A: Heat exchange algorithm 
1. Exact solution 

We would like to show that the rescaling step 

v^{t) = + {l - ()vriO) (Al) 

of the HEX algorithm is the exact solution of 

- vr), (A2) 

where K. is given by Eq. (1) and ^ by Eq. (3). We will first 
show that J^/2/C is independent of the particle velocities and 
only a function of time. This can be seen easily by considering 
the time evolution of the internal kinetic energy, which is given 
by 


heat flux. The problem with the original algorithm is that it 
exhibits a drift in the total energy whose origin remained hith¬ 
erto unclear. Eor long simulations, this energy loss becomes 
restrictive, limiting the accessible simulation time scales to a 
few nanoseconds. In our approach, we reformulated the HEX 
algorithm as a Trotter factorisation of the Liouville operator. 
Using this theoretical framework, it is straightforward to deter¬ 
mine higher-order truncation terms which are a consequence 
of the employed operator splitting. We demonstrated that the 
leading-order truncation error of the coordinates is responsible 
for the observed energy drift. 

To test the accuracy of the method, we implemented the 
eHEX algorithm in LAMMPS and ran simulations on a 
Lennard-Jones system and SPC/E water. In both cases, we 
observed at least a hundredfold reduction in the energy loss as 
compared to the HEX algorithm. With the eHEX algorithm, it 
is therefore possible to carry out constant heat flux simulations 
which are on the order of a hundred nanoseconds and based on 
fully deterministic equations of motion. 


d/C 

dt 


^ rrii (vi - vr) ■ {vi - Vp) (A3a) 

iG-y 

T (A3b) 


and therefore we can write JC{t) = JC{0) + J^t. We note that 
we can exchange the order of taking the time derivative and the 
summation, because the particle positions are fixed during this 
operation. At the same time, it is easy to see that the centre of 
mass velocity is constant in time since 
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(A4) 

(A5) 

(A6) 


In order to solve Eq. (A2) analytically, it is advantageous to 
carry out a variable transformation first. Let us consider the 
transformation Vi >—>■ Vi = Vi — vp- With the definition 
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we can express the time evolution of the new velocities as 

Vi = Xvi. (A8) 




















The solution of this equation is then given by 


It is straightforward to compute the derivative 


v,{t) = = e'u.(O). (A9) 

If we substitute the old variables back, we recover Eq. (AI) 
which proves the assertion. 


2. Splitting error 

In this section, we sketch the derivation of the leading-order 
error term of the coordinate integration arising from the oper¬ 
ator splitting in the HEX algorithm. To this end, we evaluate 
the expression 


Sri^a = (AlO) 


for the operators 
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For the first term in Eq. (AlO) we find 

[iL2, [iL2,iLi]]ri^a 

'f3,p 9 
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(Alla) 

(Allb) 


(A12a) 

(A12b) 


omitting summation bounds for /3 for readibility. In the last step 
we assumed that particles do not cross reservoir boundaries, 
in which case rji^a depends only on the velocities of particles 
within the reservoir /fc(ri)- For the second term in Eq. (AlO) 
we find 


[tLi,[iTi,*L2]]r.,„ = — V (A13) 

rUi ^ TTlj CVj^p 
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and combining the two expressions we get 


ie7f.(n) P 


dv 
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where Sij is the Kronecker delta. The final result, Eq. (20), is 
then recovered by substituting the derivative in Eq. (A 14) with 
the expression above. 
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